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Abstract 



Spatial Independent Component Analysis (ICA) decomposes the time by space functional MRI 
(fMRI) matrix into a set of 1-D basis time courses and their associated 3-D spatial maps that are 
optimized for mutual independence. When apphed to resting state fMRI (rsfMRI), ICA produces 
several spatial independent components (ICs) that seem to have biological relevance - the so-called 
resting state networks (RSNs). The ICA problem is well posed when the true data generating 
process follows a linear mixture of ICs model in terms of the identifiability of the mixing matrix. 
However, the contrast function used for promoting mutual independence in ICA is dependent on the 
finite amount of observed data and is potentially non-convex with multiple local minima. Hence, 
each run of ICA could produce potentially different IC estimates even for the same data. One 
technique to deal with this run-to-run variability of ICA was proposed by Yang et al. [2008] in 
their algorithm RAICAR which allows for the selection of only those ICs that have a high run-to- 
run reproducibility. We propose an enhancement to the original RAICAR algorithm that enables 
us to assign reproducibility p- values to each IC and allows for an objective assessment of both 
within subject and across subjects reproducibility. We call the resulting algorithm RAICAR-N 
(N stands for null hypothesis test), and we have applied it to publicly available human rsfMRI 
data (http : //www . nitre . org). Our reproducibility analyses indicated that many of the published 
RSNs in rsfMRI literature are highly reproducible. However, we found several other RSNs that are 
highly reproducible but not frequently listed in the literature. 



Notation 

• Scalars variables and functions will be denoted in a non-bold font (e.g., cr^,L,p or ^^,/). 
Vectors will be denoted in a bold font using lower case letters (e.g., y,fi,r]). Matrices will 
be denoted in bold font using upper case letters (e.g.. A, E, W). The transpose of a matrix 
A will be denoted by and its inverse will be denoted by A~^. Ip will denote the p x p 
identity matrix and will denote a vector or matrix of all zeros whose size should be clear 
from context. (^) is the number of ways of choosing L objects from N objects when order 
does not matter. 

• The jth component of vector ti will be denoted by Uj whereas the jth component of vector 
t will be denoted by tj. The element (z, j) of matrix G will be denoted by G{i^j) or Gij. 
Estimates of variables will be denoted by putting a hat on top of the variable symbol. For 
example, an estimate of s will be denoted by s. 

• If cc is a random vector with a multivariate Normal distribution with mean /x and covariance 
E then we will denote this distribution hjj\f{x \ ^i^ E). The joint density of vector s will be 
denoted by Ps(s) whereas the marginal density of si will be denoted as ps^{si). E[/(s,?7)] 
denotes the expectation of f{s^r]) with respect to both random variables s and ry. 
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1 Introduction 



Independent component analysis (ICA) [Jutten and Herault, 1991, Comon, 1994, Bell and Se- 
jnowski, 1995, Attias, 1999] models the observed data as a linear combination of a set of statisti- 
cally independent and unobservable sources. [McKeown et al., 1998] first proposed the application 
of ICA to the analysis of functional magnetic resonance imaging (fMRI) data. Subsequently, ICA 
has been applied to fMRI both as an exploratory tool for the purpose of identifying task related 
components [McKeown et al., 1998] as well as a signal clean up tool for the purpose of removing 
artifacts from the fMRI data [Tohka et al., 2008]. Recently, it has been shown that ICA applied to 
resting state fMRI (rsfMRI) in healthy subjects reveals a set of biologically meaningful spatial maps 
of independent components (ICs) that are consistent across subjects - the so called resting state 
networks (RSNs) [Beckmann et al., 2005]. Hence, there is a considerable interest in applying ICA 
to rsfMRI data in order to define the set of RSNs that characterize a particular group of human 
subjects, a disease, or a pharmacological effect. 

Several variants of the linear ICA model have been applied to fMRI data including square ICA 
(with equal number of sources and sensors) [McKeown et al., 1998], non-square ICA (with more 
sensors than sources) [Calhoun et al., 2001], and non-square ICA with additive Gaussian noise 
(noisy ICA) [Beckmann and Smith, 2004]. All of these models are well known in the ICA literature 
[Jutten and Herault, 1991, Cardoso, 1998, Comon, 1994, Attias, 1999]. Since the other ICA models 
are specializations of the noisy ICA model, we will assume a noisy ICA model henceforth. 

Remarkably, the ICA estimation problem is well posed in terms of the identifiability of the mixing 
matrix given several non- Gaussian and at most 1 Gaussian source in the overall linear mixture 
[Rao, 1969, Comon, 1994, Theis, 2004, Davies, 2004]. In the presence of more than 1 Gaussian 
source, such as in noisy ICA, the mixing matrix corresponding to the non-Gaussian part of the 
linear mixture is identifiable (upto permutation and scaling). In addition, the source distributions 
are uniquely identifiable (upto permutation and scaling) given a noisy ICA model with a particular 
Gaussian co-variance structure, for example, the isotropic diagonal co-variance. For details, see 
section 2.1.2. 

While these uniqueness results are reassuring, a number of practical difficulties prevent the reliable 
estimation of ICs on real data. These difficulties include (1) true data not describable by an ICA 
model, (2) ICA contrast function approximations, (3) multiple local minima in the ICA contrast 
function, (4) confounding Gaussian noise and (5) model order overestimation. See section 2.1.3 for 
more details. A consequence of these difficulties is that multiple ICA runs on the same data or 
different subsets of the data produce different estimates of the IC realizations. 

One technique to account for this run-to-run variability in ICA was proposed by [Himberg et al., 
2004] in their algorithm ICASSO. Using repeated runs of ICA with bootstrapped data using various 
initial conditions, ICASSO clusters ICs across ICA runs using agglomerative hierarchical clustering 
and also helps in visualizing the estimated ICs. The logic is that reliable ICs will show up in 
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almost all ICA runs and thus will form a tight cluster well separated from the rest. [Esposito et al., 
2005] proposed a technique similar to ICASSO called self-organizing group ICA (sogICA) which 
allows for clustering of ICs via hierarchical clustering in across subject ICA runs. When applied to 
multiple ICA runs across subjects, ICASSO does not restrict the IC clusters to contain only 1 IC 
from each subject per ICA run. In contrast, sogICA allows the user to select the minimum number 
of subjects for a "group representative" IC cluster containing distinct subjects. By labelling each 
ICA run as a different "subject" sogICA can also be apphed to analyze multiple ICA runs across 
subjects. 

Similar in spirit to ICASSO and sogICA, Yang et al. [2008] proposed an intuitive approach called 
RAICAR (Ranking and Averaging Independent Component Analysis by Reproducibility) for re- 
producibihty analysis of estimated ICs. The basic idea in RAICAR is to select only those ICs as 
"interesting" or "stable" which show a high run-to-run "reproducibility". RAICAR uses simple 
and automated spatial cross-correlation matrix based IC alignment which has been shown to be 
more accurate compared to ICASSO [Yang et al., 2008]. RAICAR is applicable to both within 
subject as well as across subjects reproducibility analysis. 

A few limitations of ICASSO, sogICA and RAICAR are worth noting: 

• ICASSO requires the user to select the number of IC clusters and is inapplicable without 
modification for across subjects analysis of ICA runs since the IC clusters are not restricted 
to contain only 1 IC per ICA run. 

• sogICA requires the user to select the minimum number of subjects for a "group representa- 
tive" cluster and also a cutoff on within cluster distances. 

• RAICAR uses an arbitrary threshold on the reproducibility indices selected "by eye" or set 
at an arbitrary value, such as 50% of the maximum reproducibility value. 

We propose a simple extension to RAICAR that avoids subjective user decisions and allows for 
an automatic reproducibility cutoff. The reproducibility indices calculated in RAICAR differ in 
magnitude significantly depending on whether the input to RAICAR: 

• (a) is generated using multiple ICA runs on the same data 

• (b) comes from multiple ICA runs on varying data sets (e.g. between and across subject runs) 

See Figure 1 for an illustration of this effect. Obviously, the reproducibility indices are much lower 
in case (b) since we account for both within subject and between subjects variability in estimating 
ICs. Case (b) is also of great interest from a practical point of view since we are often interested 
in making statements about a group of subjects. Hence, it is clear that a cutoff on RAICAR 
reproducibility values for the purposes of selecting the "highly reproducible" components should 
be data dependent. In this work, 

1. We propose a modification of the original RAICAR algorithm by introducing an explicit 
"null" model of no reproducibility. 
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Figure 1: Figure illustrates the variation in normalized reproducibility from RAICAR depending on 
whether the input to RAICAR is (a) Multiple ICA runs on single subject data or (b) Multiple ICA 
runs across subjects. Notice that the normalized reproducibility is much lower for across subjects 
analysis compared to within subject analysis. 
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2. We use this "null" model to automatically generate values for each IC via simulation. This 
allows for an objective cutoff specification for extracting reproducible ICs (e.g. reproducible 
at p < 0.05) within and across subjects. We call the resulting algorithm RAICAR-N (N 
stands for "null" hypothesis test). 

3. We validate RAICAR-N by applying it to publicly available human rsfMRI data. 



2 Methods 

The organization of this article revolves around the following sequence of questions which ultimately 
lead to the development of RAICAR-N: 

1. Why is a reproducibility assessment necessary in ICA analysis? In order to answer this 
question, we cover the fundamentals of ICA including identifiability issues in sections 2.1 and 
2.2. 

2. How does the original RAICAR algorithm assess reproducibility? The answer to this question 
in section 2.3 will set up the stage for RAICAR-N. 

3. How does RAICAR-N permit calculation of reproducibility p-values? In section 2.4, we 
describe the RAICAR-N "null" model and a simulation based approach for assigning values 
to ICs. 

4. How to promote diversity in group ICA runs given a limited number of subjects when using 
RAICAR-N and how to display the non-Gaussian spatial structure in estimated ICs? These 
issues are covered in section 2.5 and 2.6. 

5. How can RAICAR-N be extended for between group comparison of ICs and how does it 
compare to other approaches in the literature? This question is addressed in section 5.4. 



2.1 ICA background 

In this section, we provide a brief introduction to ICA along with a discussion of associated issues 
related to model order selection, identifiability and run-to-run variability. The noisy ICA model 
assumes that observed data y is generated as a linear combination of unobservable independent 
sources confounded with Gaussian noise: 

y ^ ^1 + As + T] (2.1) 
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In this model, 



y = p X 1 observed signal vector (2.2) 
^ — p X 1 mean vector 

A — pxq mixing matrix with p > q (more sensors than sources) and rank q 
T] = p X 1 Gaussian noise vector with density (r/ | 0, E) 
s = q X 1 vector of independent random variables (the ICs) 

with E(ss^) = D (diagonal) and E(s) = 

and with s and r] independent 

If the marginal density of the ith source Si is Psi{si) then the joint source density Ps(s) factorizes 
as Yli=iPsi{si) because of the independence assumption but is otherwise assumed to be unknown. 
Also, since the elements of s are independent their co-variance matrix D is diagonal. The set of 
variables J-" = {/x. A, £),!]} represents the unknown parameters in the noisy ICA model. Before 
discussing the identifiability of model 2.1, we briefly discuss the choice of model order or the assumed 
number of ICs q. 

2.1.1 Estimating the model order q 

Rigorous estimation of the model order q in noisy ICA is difficult as the IC densities Psi{si) are 
unknown. This means that p{y \ q^T)^ the marginal density of the observed data given the model 
order and the ICA parameters cannot be derived in closed form (by integrating out the ICs) 
without making additional assumptions on the form of IC densities. Consequently, standard model 
selection criteria such as Bayes information criterion (BIC) [Kass and Raftery, 1993] cannot be 
easily applied to the noisy ICA model to estimate q. One solution is to use a factorial mixture of 
Gaussians (MOG) joint source density model as in [Attias, 1999], and use the analytical expression 
for p{y I q^J-") in conjunction with BIC. This solution is quite general in terms of allowing for an 
arbitrary Gaussian noise co-variance E, but maximizing p{y \ q^T) with respect to T becomes 
computationally intractable using an expectation maximization (EM) algorithm for q > 13 ICs 
[Attias, 1999]. Another rigorous non-parametric approach for estimating q that is applicable to 
the noisy ICA model with isotropic diagonal Gaussian noise co- variance i.e., with SI = cr^Ip is 
the random matrix theory based sequential hypothesis testing approach of Kritchman and Nadler 
[2009]. To the best of our knowledge, these are the only 2 rigorous approaches for estimating q in 
the noisy ICA model. 

Approximate approaches for estimating q commonly used in fMRI literature (e.g., [Beckmann and 
Smith, 2004]) consist of first relaxing the isotropic diagonal noisy ICA model (with SI = cr'^Ip) into 
a probabilistic PCA (PPCA) model of [Tipping, 1999] where the source densities are assumed to 
be Gaussian i.e., where Ps{s) = J\f {s \ 0,/g). When using the PPCA model, it becomes possible 
to integrate out the Gaussian sources to get an expression for p{y \ q, J-") that can be analytically 
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maximized [Tipping, 1999]. Subsequently, methods such as BIC can be apphed to estimate q. 
Alternative approaches for estimating q in the PPCA model consist of the Bayesian model selection 
of Minka [2000], or in data-rich situations such as fMRI, even the standard technique of cross- 
validation [Hastie et al., 2009]. 

From a biological point of view, it has been argued [Cole et al., 2010] that the number of extracted 
ICs simply reflect the various equally valid views of the human functional neurobiology - smaller 
number of ICs represent a coarse view while a larger number of ICs represent a more fine grained 
view. However, it is worth noting that from a statistical point of view, over-specification of q will 
lead to over-fitting of the ICA model which might render the estimated ICs less generalizable across 
subjects. On the other hand, under-specification of q will result in incomplete IC separation. Both 
of these scenarios are undesirable. 

2.1.2 Identifiability of the noisy ICA model 

To what extent is the noisy linear ICA model identifiable? Consider a potentially different decom- 
position of the noisy ICA model 2.1: 

y = /xi + Aisi + 7/1 (2.3) 

where 

y — p X \ observed signal vector (2.4) 
/xi = p X 1 mean vector 

A\ — p X q mixing matrix with p > q (more sensors than sources) and rank q 
r]i = p X 1 Gaussian noise vector with density Af {r]i \ 0, Ei) 
si = g X 1 vector of independent random variables (the ICs) 

with E(ss^) = Di (diagonal) and E(s) = 

and with si and r/i independent 

What can be said about the equivalence between the parameterizations in 2.1 and 2.3? 

Identifiability of fi: Equating the expectations of the right hand size of 2.3 and 2.1 and noting 
that s, ?7, si, r/i have mean we get: 

Ml = M (2.5) 
Thus the mean vector /x is exactly identifiable. 
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Identifiability of A: A fundamental decomposition result states that the noisy ICA problem 
is well-posed in terms of the identifiability of the mixing matrix A upto permutation and scaling 
provided that the components of s are independent and non-Gaussian [Rao, 1969, Comon, 1994, 
Theis, 2004, Davies, 2004]. If A is a diagonal scaling matrix and P is a permutation matrix then 
the identifiability result can be stated as: 

Ai = AAP (2.6) 

where 2.3 is another decomposition of y with si containing independent and non-Gaussian com- 
ponents. In other words, the mixing matrix A is identifiable upto permutation and scaling. 

Identifiability of D and Equating the second moments of the right hand side of 2.3 and 2.1 
and noting the equality of means 2.5 and the independence of s, ij and Si, r/i we get: 

E [{y -n){y- tif] = ADA^ + S = AiD^A^ + (2.7) 

Let W he a q X p matrix and Q he a p x (p — q) orthogonal matrix such that: 

W = (A^A)-iA^ (2.8) 
Q^A^O 

Q Q ~ -^p — q 

From 2.8 and 2.7 we get: 

D + WJ:W^ = A PDiP^A^ + W^EiW^ (2.9) 
Case 1: E = a'^Ip and Ei = a^/p 

The second equation in 2.9 along with the orthogonality of Q gives = al and thus E = Ei. If 
we fix the scaling of Ai by selecting A^ = Iq then from the first equation in 2.9 we get: 

D = APDiP^A^ (2.10) 
= PDxP^A^ {PDiP^ is diagonal) 

= PDiP^ 

In other words, the noise co-variance E = cr'^Ip is uniquely determined and for a fixed scaling 
A^ Iq, the source variances D are also uniquely determined upto permutation. 

Case 2: E and Ei arbitrary positive definite matrices 

Suppose X is a square matrix and let diag(X) be the diagonal matrix obtained by setting the 
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non-diagonal elements of X to and similarly let offdiag(X) be the matrix obtained by setting 
the diagonal elements of X to 0. The noise-covariance is partially identifiable by the following 
conditions: 

Q^EQ = Q^EiQ (2.11) 
offdiag (WET^^) = offdiag {W^E^W'^) 

For a fixed scaling = /g, the sources variances D, D\ are constrained by: 

D + dmg{WT.W^) = PDiP^ + diag(TyEiTy^) (2.12) 

In general, the source variances D cannot be uniquely determined as noted in [Davies, 2004]. 

Identifiability of the distribution of s: Is the distribution of the non-Gaussian components 
of s identifiable? From 2.1 and 2.3: 

/X + A s + = /xi + Ai si + 771 (2.13) 

Substituting 2.5 and 2.6 in 2.13 we get: 

As + r] = AAPsi + r]i (2.14) 

Premultiplying both sides by W from 2.8 we get: 

s + Wt] = a P Si + Wr]i (2.15) 

Let ^^s, ^^^75 ^771 be the characteristic functions of s, r/, Si and r/i respectively. Then 

"i^sit) = E [exp (it^s)] ^^^(t) = E [exp (it^v)] (2-16) 

^s, (t) = E [exp (zt^si)] (t) = E [exp (it^m)] 

where i = \/— 1 and t is a vector of real numbers of length equal to that of the corresponding 
random vectors in 2.16. Using 2.15, we can write: 

E [exp {it^ {s + Wrj})] = E [exp {it^ {APsi + Wm})] for ah t G (2.17) 

Noting the independence of s,?7 and si,?7i: 

E [exp (it^ {s})] E [exp (zt^ {Wrj})] = E [exp {it^ {APsi})] E [exp (zt^ {Wm})] (2.18) 
^ (t) (VF^t) = (P^A^t) (VT^t) for all t G 
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Now T] and 771 are multivariate Gaussian random vectors both with mean and co- variance matrix 
SI and Hi respectively. Hence, their characteristic functions are given by [Feller, 1966, Wlodzimierz, 
1995]: 



^rj {W^t) = exp (^-^t'^W'EW^t^ for all t G 
*r,i {W^t) = exp ^-^t^WSiW^t^ for all t 



(2.19) 



Claim 2.1. A sufficient condition for identifiability upto permutation and scaling of the non- 
Gaussian distributions in s given two different parameterizations in 2. 1 and 2. 3 is: 



diagiWT.W'^) = diagiW'HiW^) 
Proof. Prom 2.20 and 2.11, we get: 

Thus from 2.19, 

{W^t) = {W^t) for all t G R'l 



(2.20) 



(2.21) 



(2.22) 



From 2.19, (W'^t) and ^ri^ {W'^t) are not equal to for any finite t, therefore, from 2.22 and 
2.18 we get: 



(f ) = ^r^^ (P^A^t) for aU t G R^^ 



(2.23) 



Note that A is a diagonal scaling matrix with entries Ai, A2, . . . , Ag on the diagonal and P is a 
permutation matrix. Thus, 



\KiiJ 



(2.24) 



where ii, ^2, ■ ■ ■ ,iq is some permutation of integers 1, 2, . . . , g. Suppose ^'s(j) is the characteristic 
function of the jth component of s and ^si(j) is the characteristic function of the jth component 
of Si. Since the components of s and Si are independent by assumption, the joint characteristic 
functions ^1/(8) and ^'(si) factorize: 



*s (*) = *s(l)(il) *s(2)(t2) ■ . . ■ • • *.(g)(tg) 



(2.25) 
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From 2.25 and 2.23 



All characteristic functions satisfy [Feller, 1966, Wlodzimierz, 1995]: 

*.W(0) = 1 (2.27) 
*si(it)(0) = 1 for all k 

Since zi, 12^ • • • ^iq is simply a permutation of integers 1, 2, . . . , there exists a j such that ij — 1. 
Then set ^2 = 0, ^3 = 0, . . . , = in 2.26. Then 2.27 and 2.26 imply: 

*s(i)(^i) = *sia)(Az,t,^.) = ^s^{j){^iti) for all ti G R (2.28) 

Select the scaling matrix as = Iq and thus A is a diagonal matrix with elements ±1 on the 
diagonal. Thus Ai = ±1 and 2.28 can be re-written as: 

*.(i)(^i) = ^s^{3){^h) for ah ti G R (2.29) 

Therefore, 

*.(i)(^i) = *.i(j)(^i) for ah G R (2.30) 
or 

^s{i){h) = *si(j)(-^i) = *-5i(j)(^i) foi" ah ti G R 

Hence the characteristic function of the 1st component of s is identical to the characteristic function 
of the (possibly sign- flipped) jth component of s\. Since characteristic functions uniquely char- 
acterize a probability distribution [Feller, 1966], the distribution of s(l) and ±si(j) is identical. 
Next, by setting ti = 0, ^3 = 0, . . . , = 0, we can find a distribution from s\ that matches the 
2nd component s{2) of s. Proceeding in a similar fashion, it is clear that the distribution of each 
component of s is uniquely identifiable upto sign flips for the choice A^ = Iq. For a general A, the 
source distributions are uniquely identifiable upto permutation and (possibly negative) scaling, as 
claimed. □ 



While the source distributions might not be uniquely identifiable for arbitrary co- variance matrices 
S, they are indeed uniquely identifiable upto permutation and scaling for the noisy ICA model 
with isotropic Gaussian noise co-variance. For more general conditions that guarantee uniqueness 
of source distributions, please see Eriksson and Koivunen [2004, 2006]. 

Corollary 2.2. IfYi — cr^Ip and Si = (Jilp, then the source distributions are uniquely identifiable 
upto sign flips for A^ = Iq. 



13 



Proof. Suppose E = a'^Ip and Si = ajlp. Then from 2.9 E = Ei and thus dmg{WT,W^) = 
diag(VFI]iVl^-^). The corohary then fohows from Claim 2.1. □ 

Corollary 2.3. If D — D\ — Iq, then the source distributions are uniquely identifiable up to sign 
flips for ^ Iq. 

Proof. U D = Di = Iq, then noting that PP^ = Iq, we get D = PD^P^ . Hence from 2.12, we 
get dmg{WT.W^) = dmg{WT.xW^). The corollary then follows from Claim 2.1. □ 



2.1.3 Why is there a run-to-run variability in estimated ICs? 

From the discussion in section 2.1.2, it is clear that for a noisy ICA model with isotropic diagonal 
additive Gaussian noise co- variance: 

1. The noisy ICA parameters F — {/x. A, Z), E} are uniquely identifiable up to permutation and 
scaling. 

2. The source distributions in s are uniquely identifiable upto permutation and scaling. 

While the above theoretical properties of ICA are reassuring, there are a number of practical 
difficulties that prevent the reliable estimation of ICs on real data: 

1. Validity of the ICA model: 

The assumption that the observed real data is generated by an ICA model is only that - 
an "assumption". If this assumption is not valid, then the uniqueness results do not hold 
anymore. 

2. Mutual information approximations: 

From an information theoretic point of view, the ICA problem is solved by minimizing a 
contrast function which is an approximation to the mutual information [Hyvarinen, 1998] 
between the ICs that depends on the finite amount of observed data. Such an approximation 
is necessary, since we do not have access to the marginal source densities Psi- Different 
approximations to mutual information will lead to different objective functions and hence 
different solutions. This is one of the reasons why different ICA algorithms often produce 
different IC estimates even for the same data. 

3. Non-convexity of ICA objective functions: 

The ICA contrast function is potentially non-convex and hence has multiple local minima. 
Since global minimization is a challenging problem by itself, most ICA algorithms will only 
converge to local minima of the ICA contrast function. The run-to-run variability of IC esti- 
mates will also depend on the number of local minima in a particular ICA contrast function. 
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4. IC estimate corruption by Gaussian noise: 

For noisy ICA, the IC realizations cannot be recovered exactly even if the true mixing matrix 
A and mean vector fi are known in 2.1. Commonly used estimators for recovering realization 
of ICs include the least squares [Beckmann and Smith, 2004] as well as the minimum mean 
square error (MMSE) [Davies, 2004]. Consider the least squares estimate s of a realization 
of s based on y : 

s = {A^A)-^A^{y - /x) = 5 + {A^A)-^A^r] (2.31) 

This means that even for known parameters, IC realization estimates s will be corrupted by 
correlated Gaussian noise. Hence using different subsets of the data under the true model 
will also lead to variability in estimated ICs. 

5. Over-fitting of the ICA model: 

Over specification of the model order leads to the problem of over-fitting in ICA. As we 
describe below, this can lead to (1) the phenomenon of IC "splitting" and (2) an increase in 
the variance of the IC estimates. 

1. IC ''splitting'' 

Suppose that the true model order or the number of non-Gaussian sources in an ICA decom- 
position of y such as 2.1 is q. Then a fundamental result in [Rao, 1969, Theorem 1] states 
that for any other ICA decomposition of y, the number of non-Gaussian sources remains the 
same while the number of Gaussian sources can change. In other words, y cannot have two 
different ICA decompositions containing different number of non-Gaussian sources. 

In view of this fact, how can a model order q ICA decomposition containing q non-Gaussian 
sources be "split" into a + 1) ICA decomposition containing {q^ 1) non-Gaussian sources 
when performing ICA estimation using an assumed model order of + 1)? As we describe 
below, the order {q + 1) ICA decomposition is only an approximation to the order q ICA 
decomposition. 

Let ai be the zth column of A in 2.1. In the presence of noise, it might be possible to 
approximate: 

aiSi^ajsj + afs^^ (2.32) 

Here: 

• aiSi is the contribution of the zth non-Gaussian source Si to the ICA model 2.1. 

• sj and sf are independent non-Gaussian random variables that are also independent 
with respect to all non-Gaussian sources sj^j ^ z in 2.1. 

• aj and a? are the basis time courses corresponding to sj and sf respectively. 
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• The time courses aj and a? look similar to each other. 

Note that if = a?, then 2.32 can be made into an equality by choosing Si = sj -\- sf . By 
replacing aiSi in 2.1 using 2.32, we arrive at an approximate model order (q+l) decomposition 
of y. In this decomposition, the component Si from a model order q decomposition appears 
to be "split" into two sub-components: sj and sf. 

2. Inflated variance of IC estimates 

Overestimation of model order will lead to over-fitting of the mixing matrix A. In other 
words, A could have several columns that are highly correlated with each other. This could 
happen as a result of IC "splitting" as discussed above. Now, for a given realization s, the 
variance of s is given by Var(s) = a'^(A'^A)~^ (for isotropic Gaussian co- variance) . An 
increase in number of columns of A and the fact that many of them are highly correlated 
implies that the variability of IC estimates Var(s) is inflated. 

In other words, running ICA multiple times on the same data or variations thereof with random 
initialization could produce different ICs. 

2.2 ICA algorithms, single subject ICA and group ICA 

In this section, we give a brief summary of how the ICA parameters are estimated in practice and 
also summarize the two most common modes of ICA application to fMRI data - single subject ICA 
(section 2.2.1) and temporal concatenation based group ICA (section 2.2.2). 

Given several independent observations y as per the noisy ICA model 2.1, most ICA algorithms 
estimate the ICA parameters T = and the realizations of s in 2 steps. We only 

consider the case with SI = cr^Ip, since as shown in section 2.1.2, the mixing matrix A and source 
distributions of s are identifiable upto permutation and scaling for this case. 

1. First, the diagonal source co- variance is arbitrarily set D — Iq. The mean vector /x is 
estimated as E(y). Then, using PC A or PPCA [Tipping, 1999], the mixing matrix A is 
estimated, upto an orthogonal rotation matrix O, to be in a signal subspace which is spanned 
by the principal eigenvectors corresponding to the largest eigenvalues of the data co- variance 
matrix E \{y — ii){y — /x)^] . The noise variance cr^ is estimated in this step as well. 

2. Next, an estimator s for the source realizations is defined using techniques such as least 
squares or MMSE. The only unknown involved in these estimates is the orthogonal rotation 
matrix O. 

3. Finally, the non-Gaussianity of the empirical density of components of s is optimized with 
respect to O using algorithms such as fixed point ICA Hyvarinen [1998, 1999]. 

For more details on noisy ICA estimation, please see [Beckmann and Smith, 2004] and for more 
details on ICA algorithms, please see [Hyvarinen et al., 2001]. 
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2.2.1 Single subject ICA 



How is ICA applied to single subject fMRI data? Suppose we are given a single subject fMRI scan 
which we rearrange as a p x n 2D matrix Y in which column i is the p x 1 observed time-course yi 
in the brain at voxel i. Observed time-courses 2/2, • • • , 2/n are considered to be n independent 
realizations of y as per the linear ICA model 2.1. Suppose S = [si, S2, . . . , Sn] is the q x n matrix 
containing the estimated source realizations at the n voxels. The jth row of S is the jth IC. In 
other words, we decompose the time by space fMRI 2D matrix into a set of basis time-courses and 
a set of q 3D IC maps using ICA. 



2.2.2 Group ICA 



How is ICA applied to data from a group of subjects in fMRI? Suppose we collect fMRI images from 
m subjects. First, we register all subjects to a common space using a registration algorithm (e.g., 
affine registration). Next, we rearrange each of the fMRI scans into m 2D matrices Yi . . . 1^, each 
of size p X n. Column j in Yi is the demeaned time-course observed at voxel location j for subject 
i. The matrices Yi . . . Ym are temporally concatenated to get a pm x n matrix Z as follows: 



(yi\ 



(2.33) 



Column i of Z is the pm x 1 vector zi which is assumed to follow a linear ICA model 2.1. 
zi,Z2,...,^n are considered to be independent realizations of the model 2.1. Suppose Sg — 
[si, S2, . . . , Sn] is a g X n matrix containing the estimated source realizations at the n voxels. The 
jth row of Sg is the jth group IC. In group ICA, the joined time-series across subjects is modeled 
using noisy hnear ICA. In practice, Yi is the PCA reduced data set for subject i. The PCA re- 
duction is either done separately for each subject using subject specific data co-variance [Calhoun 
et al., 2001] or an average data co- variance across subjects [Beckmann et al., 2005]. The aver- 
age co-variance approach requires each subject to have the same number of time points in fMRI 
scans. 



2.3 The original RAICAR algorithm 

In this section, we give a brief introduction to the RAICAR algorithm of [Yang et al., 2008]. 
Suppose we are given a data set which we decompose into nc ICs using ICA (e.g., single subject 
or group ICA). Our goal is to assess which ICs consistently show up in multiple ICA runs i.e.. 
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the reproducibility of each of these nc ICs. To that extent, we run the ICA algorithm K times. 
Suppose cc^-^^ is the n x 1 vector (e.g. spatial ICA map re-arranged into a vector) of the jth IC 
from mth ICA run. Suppose Gim is a x nc absolute spatial cross-correlation coefficient matrix 
between the ICs from runs / and m: 



GUiJ) = \corTcoef{xf\x'f^)\ (2.34) 

where |.| denotes absolute value. Gim(hj) is the absolute spatial cross-correlation coefficient be- 
tween IC i from run / and IC j from run m. The matrices Gim are then arranged as elements 
of a, K X K block-matrix G such that the Ith row and mth column of G is G(l^m) = Gim (see 
Figure 2). This block matrix G is the starting point for a RAICAR across-run component matching 
process. 

Since ICs within a particular run cannot be matched to each other, the nc x nc matrices G(/, 1)^1 — 
1 . . . K along the block-diagonal of G are set to as shown in Figure 2 with a gray color. The 
following steps are involved in a RAICAR analysis: 

1. Find the maximal element of G. Suppose this maximum occurs in matrix Gi^n at position 
(i, j). Hence component i from run / matches component j from run m. Let us label this 
matched component by MCi (the first matched component). 

2. Next, we attempt to find from each run s {s ^ I and 5 7^ m) a component that matches with 
component MCi. Suppose element {ag^j) is the maximal element in the jth column of Ggm- 
Then component is the best matching component from run s with the jth component from 
run m. 

Similarly, suppose element (i, 6^) is the maximal element in the ith row of Gis- Then compo- 
nent bg is the best matching component from run s with component i from run /. As noted 
in [Yang et al., 2008], in most cases ag — bg. However, it is possible that ag ^ bg. Hence the 
component number Cg matching MCi from run s is defined as follows: 

^ ^{cLs iiGgm{agJ)>Gig{i,bg), 
[bg if GgmiagJ) < Gig{i,bg). 

We would also like to remove component Cg of run s from further consideration during the 
matching process. To that extent, we zero out the e^th row from Ggr^r = 1 . . . K and the 
e^th column from Grs^ r — 1 . . . K . 

3. Once a matching component Cg has been found for all runs s 7^ /, m, we also zero out the zth 
row from G/^, r — 1 . . .K and the ith column from Gru r — 1 . . . K . Similarly, we zero out the 
jth column from Grm^ r — 1 . . . K and the jth row from Gmr^ r — 1 . . . K . This eliminates 
component i from run / and component j from run m from further consideration during the 
matching process. 
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4. Steps 1-3 complete the matching process for one IC component across runs. These steps 
are repeated until nc components are matched across the K runs. We label the matched 
component s as MCg which contains a set of K matching ICs one from each of the K ICA 
runs. 

Suppose matched component 5, MCg consists of the matched ICs „• . Form the 

K X K cross-correlation matrix HmCs between the matched components in MCg. The (a,6)th 
element of this matrix is simply: 

HmCs{^M = |corrcoef [xf^ ,xf^^ \ (2.36) 

The normalized reproducibility of MCg is then defined as: 



Reproducibility(MC,) = f ^—^^ j ^ HMcA^^.b) (2.37) 

^ ^ a=lb=a+l 

The double sum in 2.37 is simply the sum of the upper triangular part of HmCs excluding the 
diagonal. The normalizing factor is simply the maximum possible value of this sum. Hence 

the normalized reproducibility satisfies: Reproducibility (MCg) < 1. 

Note that our definition of normalized reproducibility is slightly different from that in Yang et al. 
[2008]. Whereas Yang et al. [2008] averages the thresholded absolute correlation coefficients, we sim- 
ply average the un-thresholded absolute correlation coefficients to compute reproducibility thereby 
avoiding the selection of a threshold on the absolute correlation coefficients. 



2.4 The RAICAR-N enhancement 

In this section, we describe how to compute reproducibility values for each matched component 
in RAICAR. Note that the RAICAR "component matching" process can be used to assess the 
reproducibility of any spatial component maps - not necessarily ICA maps. For instance, RAICAR 
can be used to assess the reproducibility of a set of PCA maps across subjects. 

In order to generate reproducibility values for the matched component maps: 

1. We need to determine the distribution of normalized reproducibility that we get from the 
RAICAR "component matching" process when the input to RAICAR represents a set of 
"non-reproducible component maps" across the K runs. 

2. In addition, we would also like to preserve the overall structure seen in the observed sets 
of spatial component maps across the K runs when generating sets of "non-reproducible 
component maps" across the K runs. 
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Gim[hj) — |corrcoef(£C^' 



(0 ^Mi 



K = number of ICA runs 

nc = number of extracted ICs in each run 

n = length of each IC 

x)l^^ = n X 1 vector of the jth IC from mth ICA run 
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diagonal nc x nc block matrices ^ 
are ignored in component matching 



maximal element of G occurs in Gim at position (i, j) 

• Hence component i from run I matches component j from run m 

• Denote this matched component by MCi 



1) 



element (ai, j) is the maximal element in the jth column of Gim 
element (i, 61) is the maximal element in the ith row of Gn 

• In this case, ai = hi. 

• Therefore we zero out the aith row from Gi^, r = 1, . . . , K 

• Similarly, we also zero out the 61th column from G^-i, r = 1, . . . , 
(zeroed out rows and columns in this case are shown in Orange) 



Add component ai from run 1 to MCi 



• element (a2, j) is the maximal element in the jth column of G2m 

• element (^, 62) is the maximal element in the ^th row of G12 

• In this case, a2 7^ 62 and G2m(Q^2,j) > Gi2{i,b2) 

• Therefore we zero out the a2th row from G2ri'r — 1, . . . , 

• Similarly, we also zero out the a2th column from Gr2i'r = 1, • 
(zeroed out rows and columns in this case are shown in Blue) 
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Add component a2 from run 2 to MCi 



After processing Grrmf l^rn and Gir., r ^ l^m we also Q 

• Zero out the zth row of Gir and ^th column of Gri for r = 1, . . . , K 

• Zero out the jth column of Grm and the jth row of Gmr for r = 1, . . . , K 
(zeroed out rows and columns in this case are shown in Green) 



Suppose matched component 1, MC\ consists of the matched ICs x^^^ , x^^^ , . . . x^^ . 
Form the K x K cross-correlation matrix HmCi between the matched compo- 
nents in MCi. The (a, 6)th element of this matrix is simply HmcA^->^) = 
|corrcoef cc^^^''^ |. The normalized reproducibility of MCi is then defined 

as: 

/ 2 \ ^ ^ 
Reproducibility (MCi) = ( ] ^ ^ Hmc, (a, h) 7 



Figure 2: Pictorial depiction of the original RAICAR algorithm [Yang et al, 2008]. The ICA 
algorithm is run K times with each run producing nc ICs. G is K xK block matrix with elements 
G{l,m) = Glm where Gim is the nc x nc absolute spatial cross-correlation matrix between ICs 
from runs / and m. The numbered green circles indicate the sequence of steps in applying RAICAR 
to a given data set. Our definition of normalized reproducibility in box 7 averages un-thresholded 
correlation coefficients thereby avoiding the selection of a correlation coefficient threshold prior to 
averaging. 



Hence for IC reproducibility assessment, we propose to use the original set of ICs across the K runs 
to generate the "non-reproducible component maps" across the K runs. 

Suppose K ICA runs are submitted to RAICAR which gives us a nc x 1 vector of observed nor- 
malized reproducibility values Reproducibility (MC^), z = l...nc - one for each IC. We propose 
to attach p-values for measuring the reproducibility of each IC in a data-driven fashion as fol- 
lows: 

1. First, we label the Knc ICs across the K runs using unique integers. In run 1, the ICs are 
labelled using integers 1, . . . , nc- In run 2, the ICs are labelled using integers (nc + 1), . . . , 2nc 
and so on. In run K, the ICs are labelled using integers (K — l)nc + 1, . . . , Knc- 

2. Our "null" hypothesis is: 

Ho : None of the ICs are reproducible (2.38) 
Hence, we can randomly label component i from run / as component d from run s 



To do this, we randomly permute the integers 1, 2, . . . ^Knc to get the permuted integers 
P(1),P(2), . . . ,p{Knc). Obviously p{i) ^ p{j) if i ^ j. 

3. The K sets "non-reproducible component runs under Hq" are constructed by assigning com- 
ponents with labels: 

• . . . ^p(nc) to run 1 under Hq. 

• p(nc + 1), . . . ,p(2nc) to run 2 under Hq 

• p ({K — l)nc + 1) , . . . , p(Knc) to run K under Hq 

4. After K runs have been generated under Hq, we subject these to a RAICAR analysis. This 
gives us nc values of normalized reproducibility, one for each matched component under Hq. 

5. Steps 1-4 are repeated R times to build up a pooled Rnc x 1 vector of normalized repro- 
ducibility Reproducibility ^^^^ under Hq. 

6. Finally, we assign a value for reproducibility to each matched IC across the K runs. The 
observed reproducibility for zth matched IC is Reproducibility(MCi) and its p-value is: 

Reproducibility^^^KMC.) = {no. of Reproducibility^,,^ ^Reproducibility + 1 

(2.39) 

7. Only those components with Reproducibilityp,^,(MCi) < PcHt are considered to be signif- 
icantly reproducible. We can use a fixed and objective value for pcrit such as 0.05. Note 
that this fixed cutoff is independent of the amount of variability in the input to RAICAR-N. 
Please see Figure 3 for a pictorial depiction of this process. 



21 



Once we calculate Reproducibility (MCi) , Reproducibility (MC2) , • ■ • , Reproducibility {MCnc)i 
how can we calculate a cutoff to separate the high reproducibihty ICs from the 
low reproducibility ICs? 



unique integer associated with 
each component in each run 



[{K -l)nc + l]...Knc 



run K 



Hq. Null hypothesis is that none of the ICs are reproducible. Hence we can 
randomly label IC i from run / as IC d from run s. This random labelling 
produces one realization of ICs under Hq . 




integers 1,2... Knc 



-^^^> permuted integers j?(l),p(2) . . .p{Knc) 



random permutation of integers 



p{l) . . . p{nc) p{nc + 1) . . . p{2nc) 



.l\Ho 



run 2 I Ho 



p{{K-l)nc + l)...p{Knc) 



run K I Hq 



reproducibility calculation for 1 realization of ICs under Hq 



• We repeat the reproducibility calculation for R realizations of ICs (e.g. 
R = 100) under Hq using the random labelling process described above. 

• This gives us a distribution of "Reproducibility" values under Hq. Each 
realization of Hq produces nc "null" reproducibility values. 

• Hence for R realizations we get Rxnc values that define the null distribu- 
tion for testing the observed reproducibility values Reproducibility(MCi). 

• Denote the vector of this "null" reproducibility values by Reproducibility^^;, 



Calculate j)- values for Reproducibility 



S ^ A 



II 



Reproducibility ^J- value for MQ = 



Observed value of Reproducibility (MC^) 

{[number of Reproducibility^^;; > Reproducibility(MCi)] + 1} 




{Rnc- 



Reproducibility 



Figure 3: Pictorial depiction of the process for generating a "null" distribution in RAICAR-N. Our 
"null" hypothesis is: "Hq: None of the ICs are reproducible. Hence, we can randomly label IC i 
from run / as IC d from run 5" . Therefore we randomly split the Knc ICs across K runs into K 
parts and run the RAICAR algorithm on each set of randomly split ICs. This gives us a set of 
"null" reproducibility values which can be used to compute values for the observed reproducibility 
of ICs in the original RAICAR run. The green circles indicate the sequence of steps for generating 
the "null" distribution after the steps in Figure 2. 
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2.5 How many subjects should be used per group ICA run in RAICAR-N? 



The input to RAICAR-N can either be single subject ICA runs or group ICA runs across a set of 
subjects. Note that the individual subject ICA runs are spatially unconstrained whereas a group 
ICA spatially constrains the group ICs across a set of subjects. Hence the number of ICs that can 
be declared as significantly reproducible at the group level are usually more than those that can 
be declared significantly reproducible at the single subject level. Hence the following question is 
relevant: 

Suppose we have a group of N subjects. We randomly select L subjects and form a single group of 
subjects. We repeat this process K times to get K groups of L subjects each of which is subjected 
to a group ICA analysis. Given the number of subjects N, how should we choose L and K? 

First, we discuss the choice of L. If L = then each of the K groups will contain the same N 
subjects and hence there will be no diversity in the K groups. We would like to control the amount 
of diversity in the K groups of L subjects. Consider any 2 subjects X and Y. The probabihty 
Pxy{L) that both X and Y appear in a set of L randomly chosen subjects from N subjects is given 
by: 

/iV-2^ 

PxYiL) = ^ (2.40) 

The expected number of times that X and Y appear together in sets of L subjects out of K 
independently drawn sets is: 

Exy{L) = K Pxy{L) (2.41) 

Ideally, we would like Exy{L) to be only a small fraction of K. Hence we impose the restric- 
tion: 



ExY (L) = KPxy{L)< amax K (2.42) 

where amax is a user defined constant such as amax = 0.05. This implies that the chosen value of 
L must satisfy: 



Pxy{L) < amax (2.43) 

In practice, we choose the largest value of L that satisfies this inequality. As shown in Figure 5, if 
iV = 23 and amax = 0.05 then the largest value of L that satisfies 2.43 is L = 5. 

The number of group ICA runs K should be as large as possible. From our experiments on real 
fMRI data we can roughly say that values of i^T > 50 give equivalent results. 
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G amma^eg, Student t and Gammapos mixture modeling 



Mixture modeling 

Figure 4: Flowchart for a group ICA based RAICAR-N analysis. The single subject data sets 
are first pre-processed and subsequently bootstrapped to create K groups, each group containing 
L distinct subjects. Each group of L subjects is submitted to a temporal concatenation group 
ICA analysis. The resulting IC maps (either raw ICs or ICs scaled by noise standard deviation) 
are subjected to a RAICAR analysis. The cross-realization cross correlation matrix (CRCM) is 
randomly permuted multiple times: G G{g,g) where ^ is a random permutation of integers 
from 1, . . . , Knc- The permuted CRCMs are subjected to a RAICAR analysis to generate a 
realization of reproducibility values under the "null" hypothesis. The computed "null" distribution 
of reproducibility values is used to assign p values to the observed reproducibility of the original 
RAICAR run. Finally, reproducible ICs are averaged using a random effects analysis and the 
resulting t-statistic images are subjected to Gamma^^e^, Student t and Gammapos mixture modeling. 
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TV = 23 subjects 




5 10 15 20 25 



L 

Figure 5: Figure shows a plot of Pxy(L) vs L for = 23 in blue. The red line shows the 
o^max = 0.05 cutoff. The largest value of L for which Pxy{L) < 0.05 is L = 5. 

2.6 How to display the estimated non-Gaussian spatial structure in ICA maps? 

The ICs have been optimized for non-Gaussianity. However, there can be many types of non- 
Gaussian distributions. It has been empirically found that the non-Gaussian distributions of ICs 
found in fMRI data have the following structure: 

1. A central Gaussian looking part and 

2. A tail that extends out on either end of the Gaussian 

It has been suggested in [Beckmann and Smith, 2004] that a Gaussian/ Gamma mixture model can 
be fitted to this distribution and the Gamma components can be thought of as representatives of 
the non-Gaussian structure. We follow a similar approach: 

1. The output of a RAICAR-N analysis is a set of spatial ICA maps (either z-transformed maps 
or raw maps) concatenated into a 4-D volume. 

2. We do a voxelwise transformation to Normality using the voxelwise empirical cumulative 
distribution function as described in [van Albada and Robinson, 2007]. 

3. Next, we submit the resulting 4-D volume to a voxelwise group analysis using ordinary least 
squares. The design matrix for group analysis depends on the question being considered. In 
our case, the design matrix was simply a single group average design. 

4. The resulting t-statistic maps are subjected to Student t, Gamma^^s and Gamma^^e^ mixture 
modeling. The logic is that if the original ICA maps are pure Gaussian (i.e., have no interest- 
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ing non-Gaussian structure) then the result of a group average analysis will be a pure Student 
t map which will be captured by a single Student t (i.e., the Gammapos and Gamma^e^ will 
be driven to class fractions). Hence the "null" hypothesis will be correctly accounted for. 

5. If the Gamma distributions have > 0.5 posterior probability at some voxels then those voxels 
are displayed in color to indicate the presence of significant non-Gaussian structure over and 
above the background Student t distribution. 

Examples of Student t, Gamma^^s and Gamma^g^ mixture model fits are shown in Figure 6. 




T score T score 

(a) (b) 



Figure 6: Examples of displaying non-Gaussian spatial structure using a Student Gamma^^s and 
Gamma^e^ mixture model. Notice how the Gamma^^e^ density is driven to near class fraction in 
the absence of significant negative non-Gaussian structure. 



3 Experiments and Results 
3.1 Human rsfMRI data 

rsfMRI data titled: Baltimore (Pekar, J . J . /Mostof sky , S.H.; n = 23 [8M/15F] ; ages: 20-40 

TR = 2.5; # slices = 47; # timepoints = 123), a part of the 1000 functional connectomes 
project, was downloaded from the Neuroimaging Informatics Tools and Resources Clearinghouse 
(NITRC) : http : //www . nitre . org/pro j ect s/f con.lOOO/. 
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Figure 7: value cutoffs for within and across single subject analysis using RAICAR-N. This figure 
illustrates the intuitive fact that within subject ICA runs are much more reproducible compared 
to across subject ICA runs. 



3.2 Preprocessing 



Data was analyzed using tools from the FMRIB software library (FSL: http://www.fmrib.ox. 
ac.uk/fsl/). Preprocessing steps included motion correction, brain extraction, spatial smoothing 
with an isotropic Gaussian kernel of 5mm FWHM and 100s high-pass temporal filtering. Spatial 
ICA was performed using a noisy ICA model as implemented in FSL MELODIC [Beckmann and 
Smith, 2004] in either single subject or multi-subject temporal concatenation mode also called 
group ICA. Please see section 2.2 for a brief summary of single subject ICA and group ICA. In 
each case, we fixed the model order of ICA at g = 40 to be consistent with the model order range 
typically extracted in rsfMRI and fMRI [Smith et al., 2009, Esposito et al., 2005]. For temporal 
concatenation based group ICA, single subject data was first affinely registered to the MNI 152 
brain and subsequently resampled to 4x4x4 resolution (MNI 4x4x4) to decrease computational 
load. 



3.3 RAICAR-N analysis with 1 ICA run per subject 



Spatial ICA was run once for each of the = 23 subjects in their native space. The resulting set 
of ICA components across subjects were transformed to MNI 4x4x4 space and were submitted to a 
RAICAR-N analysis.^ ICA components were sorted according to their reproducibility and values 

^In all RAICAR-N analyses reported in this article, we used the ^-transformed IC maps - which are basically 
the raw IC maps divided by a voxelwise estimate of noise standard deviation (named as inelodic_IC.iiii .gz in 
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were computed for each ICA component. Please see Figure 8. 




Figure 8: Single subject rsfMRI ICA runs across 23 subjects were combined using a RAICAR-N 
analysis. Figure (a) shows the observed values of normalized reproducibility (bottom) as well as 
the "null" distribution of normalized reproducibility across R — 100 simulations (top). Figure (b) 
shows the values for each IC along with the 0.05 and 0.1 cutoff lines. 

We compared the reproducible RSNs from the single subject RAICAR-N analysis to the group 
RSN maps reported in literature [Beckmann et al., 2005]. Please see Figure 9. 

To summarize, when single subject ICA runs are combined across subjects: 

• We are able to declare 4 "standard" RSNs as significantly reproducible at a value < 0.05. 

• There are 2 other "standard" RSNs that achieve a reproducibility ^- value between 0.05 and 
0.06. 

• There are 2 other "non-standard" RSNs that are of interest: one achieves a p- value of 0.0125 
and the other achieves a p-value of 0.05699. 



3.4 RAICAR-N on random sets of 5 subjects - 50 group ICA runs 

To promote diversity across the group ICA runs, as discussed in section 2.5, L = 5 subjects were 
drawn at random from the group of = 23 subjects and submitted to a temporal concatenation 
based group ICA. This process was repeated K = 50 times and the resulting set of 50 group ICA 
maps were submitted to a RAICAR-N analysis. ICA components were sorted according to their 
reproducibility and j>-values were computed for each ICA component. Please see Figure 10. 

MELODIC). It is also possible to use the raw IC maps as inputs to RAICAR-N. 
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Figure 9: The top 8 "reproducible" ICs from a RAICAR-N analysis on single subject ICA runs 
compared with standard RSN maps reported in literature [Beckmann et al., 2005]. We are able 
to declare 4 "standard" RSNs as significantly reproducible at a value < 0.05. There are 2 other 
"standard" RSNs that achieve a reproducibility value between 0.05 and 0.06 as well as 2 "non- 
standard" RSNs that achieve p-values of 0.0125 and 0.05699 respectively. We also could not find 2 
of the published RSNs in [Beckmann et al., 2005] as reproducible in single subject ICA runs. 
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Figure 10: L — b subjects were randomly drawn from the set of = 23 subjects and submitted to a 
temporal concatenation based group ICA. This process was repeated K = 50 times and the resulting 
ICA maps were submitted to a RAICAR-N analysis. Figure (a) shows the observed values of 
normalized reproducibility (bottom) as well as the "null" distribution of normalized reproducibility 
across R — 100 simulations (top). Figure (b) shows the p- values for each IC along with the 0.05 
and 0.1 cutoff lines. 

We compared the reproducible RSNs from the single subject RAICAR-N analysis to the RSN maps 
reported in literature [Beckmann et al., 2005]. Please see Figure 11. 

In summary, when 50 random 5 subject group ICA runs (from a population of 23 subjects) are 
combined using RAICAR-N: 

• We are able to declare 8 "standard" RSNs as significantly reproducible at a j)- value < 0.05. 

• There are 6 other "non-standard" RSNs that can be declared as significantly reproducible at 
a J)- value < 0.05. 



• There is 1 other "non-standard" RSN that achieves a p- value of 0.05299. 



3.5 RAICAR-N on random sets of 5 subjects - 100 group ICA runs 

To promote diversity across the group ICA runs, as discussed in section 2.5, L = 5 subjects were 
drawn at random from the group of = 23 subjects and submitted to a temporal concatenation 
based group ICA. This process was repeated K = 100 times and the resulting set of 100 group ICA 
maps were submitted to a RAICAR-N analysis. ICA components were sorted according to their 
reproducibility and p-values were computed for each ICA component. Please see Figure 12. 
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RAICAR-N Standard RSNs 




Figure 11: The top 15 "reproducible" ICs from K = 50 runs of L = 5 subject group ICA RAICAR- 
N analysis compared with standard RSN maps reported in literature [Beckmann et al., 2005]. We 
are able to declare 8 "standard" RSNs as significantly reproducible at a value of < 0.05. There 
are 6 other "non-standard" RSNs that can be declared as significantly reproducible at a value of 
< 0.05 and 1 other "non-standard" RSN that achieves a value of 0.05299. 
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Figure 12: L — b subjects were randomly drawn from the set of = 23 subjects and submitted 
to a temporal concatenation based group ICA. This process was repeated K = 100 times and the 
resulting ICA maps were submitted to a RAICAR-N analysis. Figure (a) shows the observed values 
of normalized reproducibility (bottom) as well as the "null" distribution of normalized reproducibil- 
ity across R = 100 simulations (top). Figure (b) shows the values for each IC along with the 0.05 
and 0.1 cutoff lines. 

We compared the reproducible RSNs from the single subject RAICAR-N analysis to the RSN maps 
reported in literature [Beckmann et al., 2005]. Please see Figure 13. 

In summary, when 100 random 5 subject group ICA runs (from a population of 23 subjects) are 
combined using RAICAR-N: 

• We are able to declare 8 "standard" RSNs as significantly reproducible at a value < 0.05. 

• There are 6 other "non-standard" RSNs that can be declared as significantly reproducible at 
a value < 0.05. 

• There is 1 other "non-standard" RSN that achieves a p- value of 0.05824. 
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Figure 13: The top 15 "reproducible" ICs from K = 100 runs of L = 5 subject group ICA RAICAR- 
N analysis compared with standard RSN maps reported in literature [Beckmann et al., 2005]. We 
are able to declare 8 "standard" RSNs as significantly reproducible at a value of < 0.05. There 
are 6 other "non-standard" RSNs that can be declared as significantly reproducible at a j)- value of 
< 0.05 and 1 other "non-standard" RSN that achieves a ^- value of 0.05824. 
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4 Group comparison of ICA results 



In this section, we summarize the main approaches for group analysis of ICA results which can be 
broadly classified into two categories: (1) Approaches based on a single ICA run or no ICA run 
and (2) Approaches based on multiple ICA runs. To make things concrete, suppose we have two 
groups of subjects A and B. 

4.1 Approaches based on a single group ICA run or no ICA run 

The main idea in these approaches is to use the results of a group ICA using all subjects to derive 
subject specific spatial maps for group comparison. A typical sequence of steps is as follows: 

1. The first step involves extraction of a set of template IC maps or a set of template mixing 
matrix time courses. This can be accomplished using two techniques: 

(a) Group ICA based template IC maps or time courses: 

A temporal concatenation based group ICA is run using data from all subjects in group A 
and B. This usually involves two PCA data reductions. The first reduction is based on a 
subject wise PCA decomposition [Calhoun et al., 2001] or an average PCA decomposition 
[Beckmann et al., 2005] as discussed in section 2.2.2. The next reduction is based on 
PCA reduced temporally concatenated data. Subsequently, the group ICs and the dual 
PCA reduced mixing matrix time courses are estimated using an ICA algorithm. 

(b) User supplied set of template IC maps: 

The user supplies a set of spatial maps, perhaps corresponding to an ICA decomposition 
on an independent data set. 

2. The next step either uses template IC maps or time courses. 

(a) Template time course based approach: 

First, the mixing matrix is PCA back projected and partitioned into subject specific 
sub matrices. Next, subject specific spatial maps corresponding to the group ICs are 
estimated via least-squares and a second PCA back projection is used to estimate the 
corresponding subject specific time courses. This is the approach proposed in [Calhoun 
et al., 2001], which we will refer to as the group ICA back projection approach. 

(b) Template IC based approach: 

First, spatial multiple regression using the template ICs as regressors is used against the 
original data of each subject to derive subject specific time courses corresponding to each 
template IC. Next, a second multiple regression using the subject specific time courses 
is used against the original data of each subject to derive subject specific spatial maps 
corresponding to each template IC. This approach called "dual-regression" has been 
proposed by [Beckmann et al., 2009]. A similar approach called fixed average spatial 
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ICA (FAS ICA) had also been proposed earlier in [Calhoun et al., 2004]. Both dual- 
regression and FAS ICA involve the first spatial regression stage, but dual-regression 
also includes a second temporal regression stage. 

3. Once subject specific spatial maps and time courses corresponding to group ICs have been 
determined, they are entered into a random effects analysis for group comparison. 

4.1.1 Advantages of single group ICA based approaches 

1. Much reduced computational load compared to multiple ICA based approaches. 

2. Ability to take advantage of constrained spatial IC estimation across all subjects via group 
ICA. 

Please see section 5 for discussion. 

4.2 Approaches based on multiple single subject or group ICA runs 

In these approaches results of multiple ICA runs in groups A and B are used for a between group 
analysis. A typical sequence of steps is as follows: 

1. The first step involves: 

• running a separate single subject ICA for all subjects from groups A and B (possibly 
with multiple runs per subject) or 

• running a set of group ICA runs across various sets of subjects separately, with each set 
containing subjects either from group A or group B 

2. The next step is to establish a correspondence between the ICs within and across groups. 
There are two main techniques of establishing this correspondence: 

(a) Template based methods: 

In these approaches, the user defines a template or a spatial map containing the network 
of interest. Examples of templates include a spatial map of the default mode network 
(DMN) derived from a separate ICA analysis, a spatial map from a separate PCA 
analysis, or even a binary mask defining the regions of interest. The template is then 
used to select from each run of ICA (single subject or group ICA) in each group (A 
and S), an IC that best matches the template using a predefined metric such as spatial 
correlation coefficient or goodness of fit (GOF) [Greicius et al., 2004]. 

(b) Template free methods: 

These approaches do not need a pre-defined template from the user, but instead at- 
tempt to match or cluster all ICs simultaneously within and across groups. Examples 
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of such approaches include self organizing group ICA (sogICA, [Esposito et al., 2005]) 
and RAICAR [Yang et al., 2008]. Each matched component or IC cluster includes one 
IC from each ICA run (single subject or group ICA) in each group {A and B). 

3. Finally, the selected ICs in template based methods or ICs from a selected IC cluster/matched 
component in template free methods are then entered into a random effects group analysis 
(with repeated measures for multiple single subject ICA runs) for between group comparison. 

4.2.1 Advantages of multiple ICA run approaches 

1. They account for both algorithmic and data set variability of ICA. 

2. Group comparisons happen on true ICs i.e., optimal solutions for the ICA problem. 
Please see section 5 for discussion. 

5 Discussion 

As discussed in section 2.1.2, in the noisy linear ICA model with isotropic diagonal Gaussian 
noise co-variance, for a given true model order, the mixing matrix and the source distributions 
are identifiable upto permutation and scaling. However, as pointed out in section 2.1.3, various 
factors prevent the convergence of ICA algorithms to unique IC estimates. These factors include 
ICA model not being the true data generating model, approximations to mutual information used 
in ICA algorithms, multiple local minima in ICA contrast functions, confounding Gaussian noise 
as well as variability due to model order over-estimation. A practical implication of these factors 
is that ICA algorithms converge to different IC estimates depending on how they are initialized 
and on the specific data used as input to ICA. Hence, there is a need for a rigorous assessment of 
reproducibility or generalizability of IC estimates. A set of reproducible ICs can then be used as 
ICA based characteristics of a particular group of subjects. 

We proposed an extension to the original RAICAR algorithm for reproducibility assessment of ICs 
within or across subjects. The modified algorithm called RAICAR-N builds up a "null" distribution 
of normalized reproducibility values under a random assignment of observed ICs across the K runs. 
This "null" distribution is used to compute reproducibility values for each observed matched 
component from RAICAR. An objective cutoff such as p < 0.05 can be used to detect "significantly 
reproducible" components. This avoids subjective user decisions such as selection of the number of 
clusters in ICASSO or the reproducibility cutoff in RAICAR or a cutoff on intra cluster distance 
in soglCA. 
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5.1 Results for publicly available rsfMRI data 

We applied RAICAR-N to publicly available = 23 subject rsfMRI data from http : //www . nitre . 
org/. We analyzed the data in 2 different ways: 

1. nc = 40 ICs were extracted for each of the = 23 subjects. The K = 23 single subject ICA 
runs were subjected to a RAICAR-N analysis (after registration to standard space). 

In single subject ICA based RAICAR-N analysis (see Figures 8 - 9), we are able to declare 6 
out of the 8 ICs reported in [Beckmann et al., 2005] (which used group ICA) as "reproducible" 
(4 ICs have j)- values < 0.05 and 2 ICs have values < 0.06). This is consistent with the 5 
reproducible RSNs reported in [DeLuca et al., 2005] using single subject ICA analysis. 

2. L = 5 subjects were randomly drawn from = 23 subjects to create one group of subjects 
which was subjected to a group ICA analysis in which nc = 40 components were extracted. 
This process was repeated K = 50 or 100 times and the resulting group ICA runs were 
subjected to a RAICAR-N analysis. 

In group ICA based RAICAR-N analysis (see Figures 10 - 13), we are able to declare all 8 
components reported in [Beckmann et al., 2005] as "reproducible" (at p < 0.05). Some of the 
ICs detected as "reproducible" in the group ICA based RAICAR-N on human rsfMRI data 
are not shown in [Beckmann et al., 2005] but do appear in the more recent paper [Smith et al., 
2009]. RAICAR-N results for K = 50 are almost identical to those for K = 100 suggesting 
that K = 50 runs of group ICA are sufficient for a RAICAR-N reproducibility analysis. 

5.2 Single subject ICA vs Group ICA 

Based on our results, it appears that single subject ICA maps are less reproducible compared 
to group ICA maps as illustrated in Figures 8 and 10. A single subject ICA based analysis is 
more resistant to subject specific artifacts. On the other hand, a group ICA based analysis makes 
the strong assumption that ICs are spatially identical across subjects. If this assumption is true, 
group ICA takes advantage of temporal concatenation to constrain the ICs spatially across subjects 
thereby reducing their variance. Hence, when there are no gross artifacts in individual rsfMRI data 
sets, group ICA is expected to be more sensitive for reproducible IC detection. As seen in Figures 
9 and 11, our results agree with this proposition. All ICs declared as "reproducible" in the single 
subject based RAICAR-N analysis continue to remain "reproducible" in the group ICA based 
RAICAR-N analysis. 
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5.3 How should subjects be grouped for group ICA? 

This raises the question of how the subjects should be grouped together for individual group ICA 
runs in preparation for RAICAR-N. If all N subjects are used in all group ICA analyses then there 
is no diversity in the individual group ICA runs. In this case, a RAICAR-N analysis will capture 
algorithmic variability due to non-convexity of ICA objective function but not dataset variability. 
Hence, our conclusions might not be generalizable to a different set of N subjects. 

Another option is to randomly select L subjects out of N for each group ICA run and submit the 
resulting K group ICA runs to RAICAR-N. In this case, we will account for both algorithmic and 
data set variability via a RAICAR-N analysis. In other words, we will be able to determine those 
ICs that are "reproducible" across different sets of L subjects and across multiple ICA runs. A 
key question is: How should we choose L and K7 In section 2.5, we proposed a simple method to 
determine the number of subjects L to be used in a single group ICA run out of the N subjects - 
the key idea is to form groups with enough "diversity". Multiple such group ICA runs can then be 
submitted to a RAICAR-N analysis for reproducibility assessment. Clearly, the larger the value of 
TV, the larger the value of L. Hence, increasing the number of subjects in a study will allow us to 
make conclusions that are generalizable to a larger set of L subjects. Also, conclusions generalizable 
to Li subjects are expected to hold for L2 > Li subjects but not vice versa. 

5.4 RAICAR-N for group comparisons of reproducible ICs 

In the present work, our focus was on enabling the selection of reproducible ICs for a given single 
group of subjects. However, RAICAR-N can be extended for between group analysis of reproducible 
components as well. Before we describe how to do so, it is useful to discuss other approaches for 
group analysis of RSNs described in section 5.4. Suppose we have two groups of subjects A and 
B. 

5.4.1 Discussion of single group ICA based approaches 

1. Subject specific maps corresponding to group ICA maps derived using ICA back projection 
or dual regression are not true ICs, i.e., they are not solutions to an ICA problem. 

2. These approaches do not account for either the algorithmic or the data set variability of an 
ICA decomposition. The single group ICA decomposition will contain both reproducible and 
non-reproducible ICs, but there is no systematic way to differentiate between the two. 

3. Both dual regression and ICA back projection using data derived IC templates are circular 
analyses. First, group ICA using all data is used to derive template IC maps or template 
time courses. Next least-squares based ICA back projection or dual regression using a subset 
of the same data is used to derive subject specific maps and time courses corresponding to 
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each IC. Thus model 1 (group ICA) on data V is used to learn an assumption A (template IC 
maps or template time courses) that is then used to fit model 2 (dual regression or ICA back 
projection) on a subset of the same data V. This is circular analysis [Kriegeskorte et al., 
2009, Vul and Kanwisher, 2010]. 

It is easy to avoid circular analysis in a dual regression approach via cross-validation. For 
example, one can split the groups A and B into two random parts, a "training" set and a 
"test" set. First, the "training" set can be used to derive template IC maps using group ICA. 
Next, the "training" set based template IC maps can be used as spatial regressors for dual 
regression on the "test" set. Alternatively, the template ICs for dual regression can also come 
from a separate ICA decomposition on a independent data set unrelated to groups A and 
B such as human rsfMRI data. This train/test approach cleanly avoids the circular analysis 
problem. It is not clear how to use cross-validation for an ICA back projection approach since 
template time courses cannot be assumed to remain the same across ICA decompositions. 

4. Subject specific structured noise is quite variable in terms of its spatial structure. Hence, a 
group ICA analysis cannot easily model or account for subject specific structured noise via 
group level ICs. Consequently, subject specific spatial maps in ICA back projection or dual 
regression will have a noise component that is purely driven by the amount of structured noise 
in individual subjects. On the other hand, a single subject ICA based analysis can accurately 
model subject specific structured noise via single subject ICs. 

4.2 Discussion of multiple ICA run approaches 

1. [Zuo et al., 2010] report that using different sets of template ICs in template based methods 
using spatial correlation such as [Harrison et al., 2008] can result in the selection of different 
ICs in individual ICA runs. This is not surprising since IC correspondence derived from 
template based methods does depend on the particular template used. This is similar to 
a seed based correlation analysis being dependent on the particular seed ROI used. It is 
worth noting that template free approaches such as sogICA and RAICAR do not rely on any 
template. 

2. [Cole et al., 2010] state that individual runs across subjects (or groups of subjects) can be 
quite variable in terms of the spatial structure of the estimated ICs. For example, [Cole 
et al., 2010] point out that an IC might be apparently split into two sub-components in 
some subjects but not others. The real problem is that the same model order could lead to 
over-fitting in some subjects (or groups of subjects) but not in others. Hence, the observed 
differences in a group comparison might be biased by the unknown difference in the amount 
of over-fitting across groups A and B. 

As described in 2.1.3, over- fitting can lead to the phenomenon of component "splitting" 
in ICA. This is not limited to single subject ICA but can also occur in group ICA. For 
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instance, [Zuo et al., 2010] report the "default mode" network as split into three sub networks 
using group ICA and note that component "splitting" can also reflect functional segregation 
or hierarchy within a particular IC and is not necessarily a consequence of model order 
overestimation in every case. 

Over-fitting can be correctly accounted for by a reproducibility analysis. This is because 
we expect the real and stable non-Gaussian sources to be reproducible across multiple ICA 
runs (algorithmic variability) and across different subjects or groups of subjects (data set 
variability) . 

If we want the results of a between group ICA analysis to be generalizable to an independent 
group of subjects then we must account for both the algorithmic and data variability of ICA. We 
propose to modify RAICAR-N for enabling between group comparisons of "reproducible" ICs as 
follows: 

1. Enter multiple within and across subject (or within and across sets of subjects) ICA runs 
for groups A and B into a RAICAR analysis. Perform the RAICAR component matching 
process across groups A and B. 

2. Use RAICAR-N to compute reproducibility values separately for group A and B for each 
matched component across groups A and B. 

3. Only ICs that are separately reproducible in both groups A and B and that are maximally 
similar to each other are used for between group comparisons. 

To summarize, a RAICAR-N analysis: 

%^ can be applied for "reproducible" component detection either within or across subjects in 
any component based analysis - not necessarily ICA. For instance, a set of PCA maps across 
subjects can be submitted to a RAICAR-N analysis. 

^ is simple to implement and accounts for both algorithmic and data set variability of an ICA 
decomposition. 

%^ avoids any user decisions except the final j)- value cutoff which can be objectively pre-set at 
standard values such as 0.05. 

%^ can be extended to enable comparisons of reproducible ICs between groups A and B. 



6 Conclusions 

Multiple group ICA runs using groups of subjects with enough "diversity" can be used to account 
for the run-to-run variability in ICA algorithms both due to the non-convex ICA objective function 
as well as across subjects data variability. These group ICA runs can be subjected to a RAICAR-N 
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"reproducibility" analysis. RAICAR-N enables the objective detection of "reproducible compo- 
nents" in any component based analysis of fMRI data such as ICA and can also be used for a 
between group comparison of "reproducible" ICs. 
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Figure 1: Figure illustrates the variation in normalized reproducibility from RAICAR depending on 
whether the input to RAICAR is (a) Multiple ICA runs on single subject data or (b) Multiple ICA 
runs across subjects. Notice that the normalized reproducibility is much lower for across subjects 
analysis compared to within subject analysis. 
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nc = number of extracted ICs in each run 

n = length of each IC 

x)l^^ = n X 1 vector of the jth IC from mth ICA run 



cross-realization cross-correlation matrix (CRCM) G 



Gil 






G12 




Glm 






Gn 






GiK 


























G21 






G22 










G21 






G2K 












(«2,j) 






































Gml 






Gm2 




Gmm 






Gml 






GmK 


















































Gn 






Gl2 




Glm 






Gil 






GiK 








{iM) 










































Gki 






Gk2 




GKm 






Gki 






Gkk 



diagonal nc x nc block matrices ^ 
are ignored in component matching 



maximal element of G occurs in Gim at position (i, j) 

• Hence component i from run I matches component j from run m 

• Denote this matched component by MCi 



1) 



element (ai, j) is the maximal element in the jth column of Gim 
element (i, 61) is the maximal element in the ith row of Gn 

• In this case, ai = hi. 

• Therefore we zero out the aith row from Gi^, r = 1, . . . , K 

• Similarly, we also zero out the 61th column from G^-i, r = 1, . . . , 
(zeroed out rows and columns in this case are shown in Orange) 



Add component ai from run 1 to MCi 



• element (a2, j) is the maximal element in the jth column of G2m 

• element (^, 62) is the maximal element in the ^th row of G12 

• In this case, a2 7^ 62 and G2m(Q^2,j) > Gi2{i,b2) 

• Therefore we zero out the a2th row from G2ri'r — 1, . . . , 

• Similarly, we also zero out the a2th column from Gr2i'r = 1, • 
(zeroed out rows and columns in this case are shown in Blue) 
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Add component a2 from run 2 to MCi 



After processing Grrmf l^rn and Gir., r ^ l^m we also Q 

• Zero out the zth row of Gir and ^th column of Gri for r = 1, . . . , K 

• Zero out the jth column of Grm and the jth row of Gmr for r = 1, . . . , K 
(zeroed out rows and columns in this case are shown in Green) 



Suppose matched component 1, MC\ consists of the matched ICs x^^^ , x^^^ , . . . x^^ . 
Form the K x K cross-correlation matrix HmCi between the matched compo- 
nents in MCi. The (a, 6)th element of this matrix is simply HmcA^->^) = 
|corrcoef cc^^^''^ |. The normalized reproducibility of MCi is then defined 

as: 

/ 2 \ ^ ^ 
Reproducibility (MCi) = ( ] ^ ^ Hmc, (a, h) 7 



Figure 2: Pictorial depiction of the original RAICAR algorithm Yang et al. [2008]. The ICA 
algorithm is run K times with each run producing nc ICs. G is K xK block matrix with elements 
G{l,m) = Glm where Gim is the nc x nc absolute spatial cross-correlation matrix between ICs 
from runs / and m. The numbered green circles indicate the sequence of steps in applying RAICAR 
to a given data set. Our definition of normalized reproducibility in box 7 averages un-thresholded 
correlation coefficients thereby avoiding the selection of a correlation coefficient threshold prior to 
averaging. 



Once we calculate Reproducibility (MCi) , Reproducibility (MC2) , • ■ • , Reproducibility {MCnc)i 
how can we calculate a cutoff to separate the high reproducibihty ICs from the 
low reproducibility ICs? 
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randomly label IC i from run / as IC d from run s. This random labelling 
produces one realization of ICs under Hq . 
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reproducibility calculation for 1 realization of ICs under Hq 



• We repeat the reproducibility calculation for R realizations of ICs (e.g. 
R = 100) under Hq using the random labelling process described above. 

• This gives us a distribution of "Reproducibility" values under Hq. Each 
realization of Hq produces nc "null" reproducibility values. 

• Hence for R realizations we get Rxnc values that define the null distribu- 
tion for testing the observed reproducibility values Reproducibility(MCi). 

• Denote the vector of this "null" reproducibility values by Reproducibility^^;, 
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Figure 3: Pictorial depiction of the process for generating a "null" distribution in RAICAR-N. Our 
"null" hypothesis is: "Hq: None of the ICs are reproducible. Hence, we can randomly label IC i 
from run / as IC d from run 5" . Therefore we randomly split the Knc ICs across K runs into K 
parts and run the RAICAR algorithm on each set of randomly split ICs. This gives us a set of 
"null" reproducibility values which can be used to compute values for the observed reproducibility 
of ICs in the original RAICAR run. The green circles indicate the sequence of steps for generating 
the "null" distribution after the steps in Figure 2. 
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^Resting state analysis^ 

pre-processed datasets 



Random groups of L dataset^ Random groups of L dataset^ • • • Random groups of L dataset^ 
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flCA algorithmj ' 
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G amma^eg, Student t and Gammapos mixture modeling 



Mixture modeling 

Figure 4: Flowchart for a group ICA based RAICAR-N analysis. The single subject data sets 
are first pre-processed and subsequently bootstrapped to create K groups, each group containing 
L distinct subjects. Each group of L subjects is submitted to a temporal concatenation group 
ICA analysis. The resulting IC maps (either raw ICs or ICs scaled by noise standard deviation) 
are subjected to a RAICAR analysis. The cross-realization cross correlation matrix (CRCM) is 
randomly permuted multiple times: G G{g,g) where ^ is a random permutation of integers 
from 1, . . . , Knc- The permuted CRCMs are subjected to a RAICAR analysis to generate a 
realization of reproducibility values under the "null" hypothesis. The computed "null" distribution 
of reproducibility values is used to assign p values to the observed reproducibility of the original 
RAICAR run. Finally, reproducible ICs are averaged using a random effects analysis and the 
resulting t-statistic images are subjected to Gamma^^e^, Student t and Gammapos mixture modeling. 
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L 

Figure 5: Figure shows a plot of Pxy{L) vs L for = 23 in blue. The red line shows the 
c^max — 0.05 cutoff. The largest value of L for which Pxy{L) < 0.05 is L = 5. 
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Figure 6: Examples of displaying non-Gaussian spatial structure using a Student Gammapos and 
Gamma^e^ mixture model. Notice how the Gamma^g^ density is driven to near class fraction in 
the absence of significant negative non-Gaussian structure. 
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Figure 7: value cutoffs for within and across single subject analysis using RAICAR-N. This figure 
illustrates the intuitive fact that within subject ICA runs are much more reproducible compared 
to across subject ICA runs. 




Figure 8: Single subject rsfMRI ICA runs across 23 subjects were combined using a RAICAR-N 
analysis. Figure (a) shows the observed values of normalized reproducibility (bottom) as well as 
the "null" distribution of normalized reproducibility across R = 100 simulations (top). Figure (b) 
shows the values for each IC along with the 0.05 and 0.1 cutoff lines. 
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Figure 9: The top 8 "reproducible" ICs from a RAICAR-N analysis on single subject ICA runs 
compared with standard RSN maps reported in literature Beckmann et al. [2005]. We are able to 
declare 4 "standard" RSNs as significantly reproducible at a value < 0.05. There are 2 other 
"standard" RSNs that achieve a reproducibility value between 0.05 and 0.06 as well as 2 "non- 
standard" RSNs that achieve p-values of 0.0125 and 0.05699 respectively. We also could not find 2 
of the published RSNs in Beckmann et al. [2005] as reproducible in single subject ICA runs. 
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Reproducibility cutoff 90 perc = 0.32709 




Figure 10: L = 5 subjects were randomly drawn from the set of = 23 subjects and submitted to a 
temporal concatenation based group ICA. This process was repeated K = 50 times and the resulting 
ICA maps were submitted to a RAICAR-N analysis. Figure (a) shows the observed values of 
normalized reproducibility (bottom) as well as the "null" distribution of normalized reproducibility 
across R = 100 simulations (top). Figure (b) shows the values for each IC along with the 0.05 
and 0.1 cutoff lines. 
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RAICAR-N Standard RSNs 




Figure 11: The top 15 "reproducible" ICs from K = 50 runs of L = 5 subject group ICA RAICAR- 
N analysis compared with standard RSN maps reported in literature Beckmann et al. [2005]. We 
are able to declare 8 "standard" RSNs as significantly reproducible at a value of < 0.05. There 
are 6 other "non-standard" RSNs that can be declared as significantly reproducible at a value of 
< 0.05 and 1 other "non-standard" RSN that achieves a value of 0.05299. 
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Figure 12: L = 5 subjects were randomly drawn from the set of = 23 subjects and submitted 
to a temporal concatenation based group ICA. This process was repeated K = 100 times and the 
resulting ICA maps were submitted to a RAICAR-N analysis. Figure (a) shows the observed values 
of normalized reproducibility (bottom) as well as the "null" distribution of normalized reproducibil- 
ity across R = 100 simulations (top). Figure (b) shows the values for each IC along with the 0.05 
and 0.1 cutoff lines. 
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Figure 13: The top 15 "reproducible" ICs from K = 100 runs of L = 5 subject group ICA RAICAR- 
N analysis compared with standard RSN maps reported in literature Beckmann et al. [2005]. We 
are able to declare 8 "standard" RSNs as significantly reproducible at a value of < 0.05. There 
are 6 other "non-standard" RSNs that can be declared as significantly reproducible at a j)- value of 
< 0.05 and 1 other "non-standard" RSN that achieves a ^- value of 0.05824. 
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